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ABSTRACT 



We model rapid variability of niultifrequency emission from blazars occurring across the elec- 
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fvj , tromagnetic spectrum (from radio to 7-rays). Lower energy emission is produced by the syn- 

^O ' chrotron mechanism, whereas higher energy emission is due to inverse Compton scattering of 

^•-p i the synchrotron emission. We take into account energy stratification established by particle ac- 

f-^ I celeration at shock fronts and energy losses due to synchrotron emission. We also consider the 

effect of light travel delays for the synchrotron emission that supplies the seed photons for inverse 

O , . Compton scattering. The production of a flare is caused by the collision between a relativistic 

JL I shock wave and a stationary feature in the jet (e.g., a Mach disk). The collision leads to the 

;_( ' formation of forward and reverse shocks, which confine two contiguous emission regions resulting 

C/3 . in complex profiles of simulated flares. Simulations of multifrequency flares indicate that relative 

^ ' delays between the inverse Compton flares and their synchrotron counterparts are dominated by 

energy stratification and geometry of the emitting regions, resulting in both negative and positive 

K^ , time delays depending on the frequency of observation. Light travel effects of the seed photons 

5-H ' may lead to a noticeable delay of the inverse Compton emission with respect to synchrotron 

variability if the line of sight is almost perfectly aligned with the jet. We apply the model to a 

flare in 3C 273 and derive the properties of shocked plasma responsible for the flare. We show 

that the pronounced negative time delay between the X-ray and IR light curves (X-rays peak 

after the maximum in the synchrotron emission) can be accounted for if both forward and reverse 

shocks are considered. 

Subject headings: galaxies: active — galaxies: jets — radiation mechanisms: nonthermal — 
radiative transfer — shock waves — quasars: individual (3C 273) 



1. Introduction 

In modeling emission variability from blazars, one generally assumes (e.g., Marscher 1998) that the 
spectrum is dominated by synchrotron emission at lower energies (radio through UV) and by inverse Compton 
emission at higher energies (X-rays and 7-rays). One popular model explains blazar emission as the product 



of electrons accelerated to very high energies by relativistic shocks originating somewhere near the base of 
a relativistic jet and then traveling downstream (Blandford & Konigl 1979; Marscher & Gear 1985; Hughes, 
Aller, & AUcr 1985). Emission from the accretion disk, hot corona, and molecular torus may contribute 
to the quiescent levels of observed radiation but should not affect the most rapidly variable component of 
emission, which has a time scale as short as < 1 day. The flux and level of variability of the jet is greatly 
enhanced by relativistic and light-travel effects if the angle between the jet and the line of sight is small, as 
expected in the case of blazars. 

Despite considerable observational and theoretical progress in understanding emission variability in 
recent years, there is still little agreement about a number of essential features of the problem. One of these 
is an important aspect of the inverse Compton emission: several sources of seed photons may contribute 
significantly to the observed X-ray and 7-ray flux. Among the possible sources are synchrotron photons 
produced locally or in the co-spatial regions (synchrotron self-Compton, or SSC) and external emission 
(external radiation Compton, or ERC) that can come from a variety of sources such as a molecular torus 
(Blazejowski et al. 2000), an accretion disk (Dermer & Schlickeiser 1993), and emission- line clouds (Sikora, 
Begelman, & Rees 1994). In addition, Ghisellini & Madau (1996) suggested that locally produced synchrotron 
emission may be scattered by nearby clouds surrounding the jet and enter the excited region once again where 
they can be scattered by relativistic electrons (mirror Compton, or MC). 

The location as well as the trigger of the variable high energy emission is currently an open question. 
Propagating shock waves can produce long-term (months to years) outbursts (Marscher & Gear 1985; Hughes, 
Aller, & Aller 1985). At high frequencies, interaction with a turbulent ambient medium in the jet can generate 
more rapid fluctuations (Marscher & Gear 1985; Travis 1992). The onset of a flare could also be triggered 
by the collision between a relativistic shock wave traveling down the jet and a more slowly moving (Spada et 
al. 2001) or stationary compression, e.g., the core seen in very long baseline interferometric (VLBI) images 
of blazars. Here we consider this last possibility, that a moving shock collides with a stationary formation 
associated with the core of the jet, which results in rapid variability on the time scale of several days. 

In order to progress with such models, the nature of the radio core needs to be defined. Synchrotron 
self- absorption determines the location of what appears to be the core at lower radio frequencies, manifested 
by an upstream shift in position as one increases the frequency of observation. However, at sufficiently 
high radio frequencies (> 40 GHz), we expect this trend to cease so that images reveal the "true" core that 
corresponds to a stationary physical structure in the jet. 

Cross-correlation analysis of long term variability of such blazars as PKS 1510—089, 3C 273, and 3C 279 
reveals strong correlations between high (X-rays) and low (optical, IR, and radio) energy emission, suggesting 
that both synchrotron and inverse Compton emission arise from the same location in the jets (McHardy et 
al. 1999; Marscher et al. 2004). Moreover, a comparison of multi-epoch VLBI images and long term light 
curves of these blazars suggests that the flares are associated with activity in or near the radio cores. Bright 
flares usually correspond to an increase in brightness of the radio core followed by the appearance of a 
superluminal feature emerging and traveling downstream from the core. This analysis supports the notion 
that the flares may be generated by a collision between a moving shock wave and a stationary structure. 

Theoretical considerations (Courant & Friedrichs 1948) as well as a number of numerical simulations 
of relativistic jets (Dubai & Pantano 1993; Gomez et al. 1997) have shown that stationary structures may 
appear in jets as a result of interaction between the ambient medium and the jet if there is local pressure 
imbalance between the two. This interaction leads to the formation of a system of standing oblique shock 
waves, perhaps terminating in a strong shock (Mach disk) that lies perpendicular to the jet flow. A relativistic 



flow of plasma crosses the Mack disk, which accelerates the plasma electrons to highly relativistic energies. 
Synchrotron radiation produced by these particles then gives rise to the emission seen as the radio core in 
VLBI images. 

In this paper, we present a theoretical model for rapid emission variability from blazars based on 
the scenarios sketched above. We then discuss the properties of the model that can help constrain the 
physical properties of the emission regions. We apply the model to a flare from the quasar 3C 273 observed 
simultaneously in the infrared K-band and in X-rays with RXTE in March 1999. We fit the light curves 
and the evolution of the X-ray spectral index during the flare and determine the physical parameters of the 
medium producing the emission. 



2. Previous Observational and Theoretical Work 

2.1. Multifrequency Light Curves and VLBI observations 

Since 1996, Marscher et al. (2004) have been conducting multifrequency monitoring of the quasars 
PKS 1510-089 and 3C 279 with RXTE in the 2.4-20 keV X-ray band as weh as in several other spectral 
bands, including radio and optical. The quasar 3C 273, which exhibits some blazar characteristics, was also 
observed simultaneously in the X-rays and the infrared K-band during concentrated campaigns in January 
1997 (McHardy et al. 1999) and March 1999. Cross-correlation of the light curves shows that there is 
a correspondence in variability between different bands. This is particularly pronounced in two flares in 
PKS 1510-089 observed in 1997. In addition, violent activity of 3C 279 in 1998 and 2001 in X-rays is wefl 
matched by contemporaneous variations in the optical. There is also a general correspondence between the 
ejection of supcrluminal components from the core and the onset of the flares. 

Stressing parallels between AGNs and microquasars in our galaxy, Marscher et al. (2002) conducted a 
study of long-term X-ray light curves and VLBI maps of 3C 120. They determined a minimum distance of 
0.3 pc between the X-ray source — identified with the accretion disk for this particular galaxy — and the core 
of the radio jet. The radio core is likely to be even farther downstream for higher luminosity AGNs that 
harbor more massive black holes in their nuclei. 

These observations suggest a model of blazar variability in which both synchrotron (optical, IR, radio) 
and inverse Compton emission (X-rays) originate from the same region or contiguous regions, which are 
located in the relativistic jet at a site near the radio core seen in VLBI images. It is natural to assume that 
the onset of a flare is triggered by the collision between a moving shock wave and the stationary structure 
in the jet (e.g., Mach disk) responsible for the radio core. The excited plasma resulting from the collision 
moves down the jet and is later seen as a supcrluminal feature. We develop such a model in Sections 3-6 
below. 



2.2. Acceleration of particles by relativistic shock w^aves 

Radio to infrared flares in blazars are usually explained by appealing to the shock in jet model (Blandford 
& Konigl 1979; Marscher & Gear 1985; Hughes, Aller, & Aller 1985), which postulates that relativistic shock 
waves in the jet cause particle acceleration and nonthermal emission, appearing as supcrluminal knots in 
VLBI maps. Relativistic shocks in a jet may result from intermittent variability in the physical properties 
of the portion of the central engine that is connected with the formation of the jet. The details of the 



jet production and coUiniation have been associated with differential rotation in the inner accretion disk, 
which creates a magnetic coil that expels accreting material (Meier, Koide, & Uchida 2001). It is likely 
that the quasi-steady process that leads to jet production is sometimes punctuated by dramatic events 
due to instabilities occurring in the accretion disk surrounding the black hole. One of these events may 
momentarily saturate the jet with extremely energetic plasma with much higher pressure than the steady 
jet plasma downstream. The high energy plasma will sweep up slower moving material, forming a shell of 
enhanced density that is separated from the steady jet in front of it by a shock wave. This is usually referred 
to as a "forward shock" since the shock front advances ahead of the shocked material. On the upstream 
side, the shell is separated from the high energy plasma that pushes it by a contact discontinuity. Inside 
the energetic plasma, the information about its interaction with the slower flow ahead is transmitted by a 
compression wave that can steepen into a reverse shock. The exact outcome depends on the details of the 
jet production and the instability that injects energetic plasma into the steady jet (BickncU & Wagner 2002; 
Marti, Miillcr, & Ibaiiez 1994). 

A Fermi-type acceleration mechanism is typically invoked to explain how particles gain the high energies 
required by observations of the synchrotron and inverse Compton emission (see, e.g., Krolik 1999). The effects 
of relativistic beaming restrict the conditions under which the repeated cycles of acceleration may proceed, 
but the resulting energy spectrum of the accelerated particles is similar to that observed in the non- relativistic 
case. Achterberg et al. (2001) performed stochastic simulations of the acceleration process at ultrarelativistic 
shock fronts and found that a power law spectrum N{E) ex E~'' is produced with a nearly universal value of 
the slope, s « 2.2-2.3. The maximum energy attainable in the acceleration process is limited by the particle 
diffusion that eventually puts it outside of the restricted range of conditions needed for the acceleration 
cycles to proceed. Moreover, extremely high energy particles produced by acceleration should suffer from 
severe radiative energy losses through synchrotron and sometimes inverse Compton emission. Radiative 
losses can steepen the energy spectrum of relativistic particles above a certain break energy and also limit 
the maximum value of the energy in the spectrum. 



2.3. Stationary shocks in Gas Dynamical Simulations of Jets 

A stationary shock wave system in a relativistic jet can be caused by pressure imbalance between the 
jet and the surrounding ambient medium. Oscillations of the jet width occur as the overpressured medium 
in the jet expands until its pressure falls below the ambient gas pressure. This underpressured plasma then 
contracts under the influence of the external medium until high pressure is restored. This sets up a stationary 
system of oblique shocks and a Mach disk if the pressure mismatch is substantial (see Courant & Friedrichs 
1948). 

Daly & Marscher (1988) used the method of characteristics in two dimensions to determine that standing 
shocks appear in relativistic jets in a similar fashion to those in non-relativistic flows. They speculated that 
the radio core is produced by plasma passing through an oblique shock/Mach disk system in the jet. In 
this scenario, stationary hot spots seen downstream of the core (Jorstad et al. 2001) may be explained by 
the oscillatory behavior of the jet producing multiple sets of stationary structures. Dubai & Pantano (1993) 
also employed the method of characteristics in order to study the steady-state structure of relativistic jets 
with toroidal magnetic field. They found that for constant external pressure the jet has strong periodic 
features with large pressure and velocity changes. These authors noted that the features could be associated 
with a locally strong toroidal magnetic field. In their simulations, the strength of the features downstream 
decreased substantially when the radial pressure gradient of the external medium was introduced. Using a 



numerical relativistic fluid dynamics code, Komissarov & Falle (1997) modeled the appearance of jets on 
parsec scales by including moving shocks and a stationary reconfinement shock. They identified moving 
shocks in their simulations as the observed superluminal knots and the reconfinement shock as the brightest 
stationary knot, i.e., the radio core. 

Gomez et al. (1997) conducted hydrodynamical simulations of the generation, evolution, and radio 
emission of superluminal components in axially symmetric relativistic jets. Their simulations showed that 
an increase in speed of the flow at the jet base results in a relativistic shock wave propagating down the 
jet. Prior to the introduction of velocity perturbations, their calculations produced a steady jet that showed 
a periodic system of stationary oblique shocks when the jet was initially overpressured with respect to the 
external medium. Gomez et al. (1997) notice that when the moving perturbation passes through a standing 
shock the latter is dragged downstream for some distance before returning to its initial position as the steady 
jet is reestablished. Three-dimensional hydrodynamic simulations of relativistic jets (Aloy et al. 2003) show 
more complex behavior when precession of the jet is introduced. 

According to the simulations by Gomez et al. (1997) the Mach disk is expected to have a small size 
compared to the cross-section of the jet. Most of the plasma in the jet interacts with the oblique shocks, 
which decrease the bulk Lorentz factor by a moderate amount, after which the fiow reaccelerates. Even 
the plasma that passes through the Mach disk radiates only a small fraction of its internal energy before 
it reaccelerates into the lower density region downstream and once again reaches relativistic speed. It can 
do so as long as the slope of the electron energy spectrum is steeper than 2 or if the protons carry most 
of the kinetic energy. Otherwise, a large fraction of the total energy will be dissipated through radiation 
by the highest energy electrons. Therefore, the radio core in most sources is not expected to be nearly as 
intrinsically luminous as the hot spots in the external radio lobes where most of the relativistic energy of 
the jet is deposited. In some sources, notably TeV blazars, most of the energy in the jet indeed seems to be 
dissipated on parsec scales, which is an expected outcome of our model if the distribution of electrons is fiat 
with slope s < 2 and the protons do not carry the bulk of the total energy out to parsec scales. 

The association of the radio core with the underlying standing structure in the jet is complex. It cannot 
be equated with the plasma slowed down and compressed by the Mach disk since this is very weakly beamed. 
VLBI observations of single radio cores suggest that the emission from radio cores is beamed in the direction 
along the jets. The sources with jets in the plane of the sky do not possess strong radio cores, which again 
suggests relativistic beaming. The contributions from oblique shocks and also from the plasma shocked by 
the Mach disk after it reaccelerates to relativistic speeds must be taken into account. 



2.4. Progress in Modeling 

In this section we review recent progress in theoretical understanding of multifrequency variability of 
blazars. We pay special attention to the basic assumptions made in the simulations about the geometry, 
acceleration, and internal structure of the emitting regions. 

In simulating the spectra of 7-ray emitting blazars, Mastichiadis & Kirk (1995) (see also Mastichiadis & 
Kirk 1997) assumed that the emission originates in a spherical, homogeneous region into which relativistic 
electrons with a power law energy distribution are injected. They approximated that the injection has a 
fixed form without identifying a specific acceleration mechanism but interpreted the overall picture as one 
in which a shock front rapidly accelerates the electrons. The light-travel time of the photons across the 
emitting region was neglected. Changes in maximum energy of electrons Emax, magnetic field _B, or the 



normalization factor in the distribution of electrons Nq, were considered by these authors as the factors that 
cause the flares. 

Dermer et al. (1997) modeled flares by computing the time dependent spectral flux of a relativistic, 
homogeneous, spherical blob of plasma with instantaneous injection of relativistic electrons. The flares 
considered were restricted to ones in which the energy loss of relativistic electrons is dominated by synchrotron 
or external Compton emission. The external radiation that drives the inverse Compton scattering in the 
ERC model was assumed to be isotropic, homogeneous, and constant, similar to what is expected when 
the emitting blob is inside the broad emission line region (BLR). A (5- function approximation for Compton 
scattering and synchrotron emission was used in computing the observed fluxes. Energy losses due to 
adiabatic expansion were neglected in this work. Dermer (1998) presented the results of the modeling of 
multifrequency flares based on the approach presented in Dermer et al. (1997), stressing the correlations 
between optical and X-ray emission. 

Coppi & Aharonian (1999) presented sample SSC calculations based on a one-zone homogeneous model 
developed by Coppi (1992) that mostly concentrated on an improved treatment of the microphysics at 
the expense of more relevant but poorly understood dynamical and structural considerations. A delay in 
Compton flares observed in the simulations was explained by Compton radiation at relatively low frequencies 
being emitted by lower energy electrons with longer lifetimes compared to higher energy electrons that 
produce high frequency synchrotron emission. They also note that in the SSC model the target photon 
intensity in the source can only change on a time scale defined by the light travel time through the source 
or longer. 

Chiaberge & Ghisellini (1999) pointed out that if the electron distribution changes on a time scale 
shorter than the light crossing time through the emission region, then even a homogeneous source, i.e., 
one with constant density of relativistic electrons throughout, will resemble an inhomogeneous one. They 
stressed that considering only the time evolution of the energy distribution of the emitting electrons is not 
enough for correctly modeling the flares. Due to light travel effects, the observer will see radiation produced 
in separate parts of the source at different times and, therefore, characterized by particle distributions having 
various ages. Both simultaneous injection of electrons through the entire source and energization occurring 
at a shock front traveling through a section of the jet were considered. The emitting region was divided into 
parallel slices to account for light travel effects. In the case of the shock injection, each slice was further 
divided into separate cells occupied by particles of different age. The angle between the jet axis and the 
direction toward the observer was assumed to he 9 ^ 1/r, where F is the bulk speed of the shocked plasma. 
The authors noted that the inverse Compton emission calculation requires consideration of light travel effects 
as well, since the seed photons in the SSC model may travel across the source on the same time scale as the 
scale of variability. This could lead to a time delay of the Compton emission with respect to the synchrotron. 
They did not, however, include the effect in their calculations. The authors also stressed that the symmetric 
shapes of the light curves strongly constrain the injection and cooling time scales. 

Kataoka et al. (2000) computed time dependent SSC models to study the spectral evolution and variabil- 
ity patterns of blazars based on a one-zone homogeneous approximation with electrons injected uniformly 
throughout the emission region. They quantitatively compared the model with observed X-ray flares to 
derive the physical parameters of the emitting region. However, only the X-ray bands were used because 
the source was not observed in other spectral regions during the flare. Similar to the work by Chiaberge & 
Ghisellini (1999), the authors included the effects introduced by radiation travel times of photons originating 
in different parts of the emission region, but neglected the light travel delays of the seed radiation for the 
inverse Compton scattering. 



An illustration of the importance of light travel delays in estimating the observed variability of emission 
from blazars was demonstrated in Salvati, Spada, & Pacini (1998), who suggested a model for explaining 
an extremely rapid burst observed from Mrk 421. In their model, the observed radiation is produced as the 
result of prolonged interaction between a disturbance traveling down the jet and conical shocks. Due to the 
specific conical geometry of the shocks chosen by the authors, the photons that are emitted in the observer's 
direction arrive at the observer at essentially the same time, producing a rapid, extremely intense burst of 
radiation. 

Some models (see Mcszaros 2002) interpret gamma-ray bursts (GRBs) as shells of plasma ejected from a 
center site, colliding as they move outward. A system of forward and reverse shocks is established as a result 
of the collision. These shocks accelerate plasma particles, which in turn radiate their energy and produce a 
multifrequency flare. Spada et al. (2001) applied this "internal shock" scenario to blazar variability. They 
attributed the formation of shells to intermittent variability of the central engine. Plasma characteristics 
were treated not as parameters but rather as the result of the jet dynamics. Despite appealing to forward 
and reverse shock waves, the authors assumed that the emitting zone is homogeneous with the same energy 
distribution throughout the entire emitting volume. Light travel delays across the emitting region were 
neglected in their study. 

In modeling broad-band spectra of blazars that include SSC emission as well as ERG from both the BLR 
and the dusty torus surrounding AGN, Blazejowski et al. (2000) found that the BLR may be less significant 
as the source of the seed photons than the radiation from hot dust in the torus if the emission region in the 
jet is located at a distance greater than ~ 0.1 pc from the central engine. Reverberation mapping techniques 
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(Peterson 1993) indicate that the average radius of the BLRs is ~ O.IL42 pc, where L42 is UV luminosity 
normalized by 10''^ erg s^^ A . Therefore, if the emitting shell in the jet is outside the BLR, then the 
external radiation from emission line clouds will be Doppler boosted by a lesser amount compared to the 
boosting of emission from the torus. The results of the modeling (Blazejowski et al. 2000) showed that the 
ERG light curve may still be dominated by the contribution from the BLR at the highest observed 7-ray 
energies simply due to the higher frequencies of the BLR photons compared with the IR photons from the 
dust. Further modeling by Sikora et al. (2001) gave strong support to the SSG mechanism as the primary 
emission process in the X-ray band. 



3. Gas Dynamics of Shock Collision 

In our model, calculations of the production and transfer of the synchrotron and inverse Gompton 
emission are performed in the rest frame of the emitting plasma. The input parameters that describe the 
geometry of the source and the state of the plasma are defined in the plasma rest frame as well. If the bulk 
speed of the emitting plasma is given, the observed quantities such as synchrotron and inverse Gompton flux 
as a function of frequency and time can be calculated provided that the angle between the jet axis and the 
line of sight toward the observer is known. Since variability of blazars seems to be associated with the radio 
core, we assume that the collision between a moving shock and a stationary shock is the primary source of 
multifrequency flares and outbursts. In this section, we describe how the properties of the shocked plasma 
relate to the conditions in the quiescent jet flow and the colliding shocks. 

We consider a collision between a shock moving down the jet at speed v'g and a stationary Mach disk 
(see Fig. 1). All the primed quantities are evaluated in the AGN rest frame, whereas all the plasma rest 
frame quantities are unprimcd. Prior to collision, the moving shock wave approaches stationary Mach disk. 



which results in gradual shrinking of zone (3) and eventually a head-on collision between the moving shock 
and the Mach disk. The reverse and forward shocks (separated by a contact discontinuity) that result from 
the collision arc shown in the space-time diagram in Fig. 2. Both shocks propagate downstream in the frame 
of the AGN, but move in opposite directions in the shocked plasma rest frame. Shocked plasma occupies 
zones (r) and (/). Following the collision that occurs at time i' = 0, the forward shock expands into the 
dense medium accumulated behind the Mach disk [zone (4)] until it reaches the rarefaction region [zone (5)] 
at time t'^f. Similarly, the reverse shock expands into zone (2), which consists of the material swept up by 
the moving shock prior to collision. We assume that the acceleration in the reverse zone effectively stops 
once the reverse shock reaches the contact discontinuity between (1) and (2) at time t'^^. It should be noted 
that the forward and reverse shocks interacting with the contact discontinuities or the rarefaction behind the 
Mach disk may set up secondary shocks and/or waves that could complicate the system further (see Gomez 
et al. 1997); we ignore these affects for the sake of simplicity. 

The properties of the shocked plasma can be found by using the shock jump conditions for a relativistic 
gas (see Anile 1989), which are based on the conservation of matter, momentum, and energy across the 
shock front. In this section we set the speed of light c to unity for convenience. For ID flow, the conserved 
quantities reduce to 

pT{v — v^) = const, 
T'^{v — v^){e+p) +pvj: = const, (1) 

r {v — v-s)(e + p)v + p = const. 

These relations connect proper density p, pressure p, internal energy e, and plasma speed v in units of c (the 
corresponding Lorentz factor F — l/\/l — v-^) on both sides of the shock front, which moves at speed vs. 
They can be applied in any reference frame. If conditions upstream of the shock as well as the speed of the 
shock vs are known, the shocked plasma parameters can be calculated provided that the equation of state 
is specified. We assume a fully relativistic plasma, in which case the equation of state is p = e/3. 

Let us assume that plasma properties (/03,P3,f3) in zone (3) are known. In order to make the problem 
tractable, we approximate that velocity, pressure, and density are constant within each zone. Since the 
speeds of both the moving shock and the Mach disk are known (v'j^j ~ for the Mach disk) , we can therefore 
determine all the plasma parameters in zones (2) and (4) as described below. However, the speeds of 
the forward and reverse shocks are unknown, so the conditions in zones (r) and (/) must be determined 
simultaneously by appealing to the constancy of the plasma speed v' and pressure p across the contact 
discontinuity between (r) and (/). Applying shock jump formulas to both zones in the plasma rest frame 
gives six equations for six unknowns: p^, Pf, Vr, Vf, p, and v' . After some substitutions we find the following 
equation: 

1 — 2iVrV2 1 — 3u/W4 



P2(1+W,.W2) Piil+VfViY 



(2) 



where Vf = {-s/va + 3 — V4)/i and Vr — {\/v2 + 3 — W2)/3 are the forward and reverse shock speeds in the 
plasma rest frame. In this reference frame the shock speeds cannot exceed l/v3 but they are also greater 
than 1/3 for a fully relativistic plasma. Eq. (2) can be solved for v' by applying relativistic transformation 
formulas, 

Vi = z TT and V2 = ^^, (3) 

1 - v'pV'^ 1 - v'^v'p 

that connect the AGN rest frame speeds (primed quantities) and the plasma rest frame speeds (unprimed). 
Since v'^ < v' < V2 in the AGN rest frame, speeds V2, V4, Vf, v^ are defined to be positive in the above 



formulas for convenience despite the fact that the shocks move in opposite directions in the plasma rest 
frame. 

Although Eq. (2) depends on p2 and p4, these quantities scale as ps for a fully relativistic plasma and 
therefore the plasma speed v' only depends on the quiescent flow speed Wg and the moving shock speed 
v'g. In general, Eq. (2) must be solved numerically, but an important analytical solution can be obtained 
if p2 — Pi- This is equivalent to stipulating that in the rest frame of the quiescent flow [zone (3)] the 
speeds of the moving shock and the Mach disk are equal. This condition is satisfied if v'^ = 2f3/(l + w'g) or 
r^ = r'3(I + v' ^). The resulting shocked plasma speed is v' = Wg, which also gives p,. = pf and Vr = f/. 
Indeed, if the colliding shock speeds are equal in the rest frame of the quiescent flow then the shocked plasma 
will remain stationary in this frame while the speeds of the reflected shocks will be equal. 

Fig. 3 presents the shocked plasma Lorentz factor as a function of the moving shock Lorentz factor (top 
panel) for a constant value Fg = 4 and pg = 1, pg = 1. We have used the equation of state for a fully 
relativistic plasma in all zones. The minimum possible shock speed of Vg — l/\/3 in the quiescent plasma 
rest frame, combined with the bulk Lorentz factor of the quiescent plasma, give the lower limit on possible 
values of the moving shock Lorentz factor, F^ ^j„ « 7.6. The condition p2 = Pi is satisfied when F'^ « 31 for 
the chosen value F' = 4. 

Using the analytical solution v' = Ug discussed in the previous paragraph allows one to find an approx- 
imation, 

F' 

r; "-^, (4) 

that gives an exact solution for the p2 = Pi case and does not deviate much from the exact solution when 
Fg ~ Fg (1 -l-Wg ) and F^ :^ 1. Approximation (4) is shown in the top panel of Fig. 3 (dotted fine). It should 
be noted that extremely fast relativistic shocks are required to achieve moderate shocked plasma speeds v' 
e.g., F^ = 31 produces only F' = 4. 

Reverse and forward shock speeds and plasma densities in the shocked plasma rest frame are shown in 
the middle and bottom panels of Fig. 3. An interesting result is that the value of the reverse and forward 
densities pr and p/ are similar for most values of F^ and Fg (see bottom panel of Fig. 3). This is to be 
expected when p2 ^ Pi- To understand why the densities are still close even when p2 is much different from 
P4, consider that the final densities pr and pf result from two consecutive amplifications from the starting 
value of p3. The two amplifications that result in the final density pf occur at the Mach disk and the forward 
shock, and at the moving shock and the reverse shock for pr- For example, F'^ = 9 is well below the critical 
value of 31 that gives p,- = pf (Fig. 2). Thus, the density amplification caused by the moving shock is 
small compared with that caused by the Mach disk. Since the moving shock is relatively weak, its collision 
with the Mach disk will slow down the plasma to a mildly relativistic speed (see Fig. 2). Therefore, the 
amplification by the forward shock will be small whereas that of the reverse shock will be relatively large 
because it encounters fast plasma from zone (2). For this reason, the forward and reverse densities resulting 
from these double amplifications are similar. The same reasoning applies to the case when the moving shock 
is stronger than the Mach disk, resulting in faster moving plasma. The amplification of the transverse (to 
the flow) component of dynamically insignificant magnetic field should exhibit similar behavior. 

We note that the above is strictly true only if quiescent-flow plasma in zone (3) is homogeneous, without 
any gradients along the jet axis. Since the jet may be expected to converge somewhat at the position of the 
Mach disk, pg must be larger at the Mach disk. If properly taken into account, this could result in higher 
values of the forward density pf and magnetic field strength Bf. 
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4. Structure of the emission region 

It is clear that once the profiles of the individual flares are resolved, a more detailed model of the emission 
region that takes into account the excitation structure, geometry of the region, and the finite speed of light 
propagation must be used to reproduce observed flares. It is usually assumed that the relativistic electrons 
are injected uniformly in space within the source. The mechanism that can provide this type of injection is 
often left unspecified, and it is indeed difficult to conceive of a physical process that can change conditions 
simultaneously across the emission region. We therefore adopt a scenario in which changes are transmitted 
by an excitation front, which we assume to be a shock. Because of highly relativistic motion of the source 
in the observer's frame, one needs to take into account the light travel time and relativistic effects that may 
dramatically influence time scales, frequencies, and flux levels of the observed radiation. Furthermore, in 
the SSC model another type of light-travel time effect must be considered. The inverse Compton scattering 
at a particular location within the confines of the excited region is likely to involve seed photons produced 
throughout the region. Many of these seed photons travel across a significant fraction of the excited region. 
These travel times may be comparable to the characteristic time scale of structural changes of the source 
excited by a relativistic shock wave moving at almost the speed of light. 



4.1. Excitation structure and geometry 

In this section we discuss the geometry of the source, the excitation structure introduced by a passing 
shock wave, and the effect of light travel delays on the internal structure as viewed from different locations, 
for example, as seen by a distant observer. In our model, we approximate that the excitation region has a 
cylindrical shape with the axis of the cylinder coinciding with the axis of the jet. Its transverse extent is 
determined by the size of the Mach disk, which could be much smaller than the cross-section of the jet itself. 
The region is filled uniformly with relativistic electrons gyrating in a mostly turbulent magnetic field. 

A snapshot of the excitation structure after the collision is shown in Fig. 4. The electrons are accelerated 
in a very narrow region at the shock wave fronts, which are assumed to be much smaller than the dimensions 
of the source. As soon as electrons leave the acceleration region, their energy steadily decreases due to 
radiative energy losses. Although the density of electrons does not change, the maximum energy of electrons 
declines, establishing a gradient in the electron distribution across the source. The farther away from the 
shock front a given location is, the lower is the maximum energy in the distribution. The excitation structure 
follows the progression of shocks through the forward and reverse regions, which are assumed to be moving 
at constant speeds. 

Because of intrinsic similarity between the forward and reverse regions, here and in the following sections 
we mostly consider the formulas for the forward emission region. We drop the subscript / unless both 
regions must be considered. Analogous expressions for the reverse region can be obtained straightforwardly 
by inverting the z-axis. The evolution of the distribution of electrons is found by solving the equation 

with the initial distribution N{j) = No^~'^ and 7m,m < 7 < Jmax- We assume that there are no sources 
or sinks of electrons in the region between the acceleration zone and outer edge of the volume of excited 
electrons. We also neglect adiabatic losses since the expansion of the emitting volume is insignificant during 
relatively fast flares caused by collisions with the Mach disk. The solution depends on time t and position 
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along the jet axis z. The energy loss function 7 depends on 7 and the average magnetic field strength B if 
synchrotron losses provide the dominant contribution: 7 = — 7^/ii, where t\ = 7.73 x 10* s/(_B/l G)^. If 
inverse Compton losses are significant, then 7 will also depend on the properties of the seed photon field. In 
the case of the SSC model, the inclusion of Compton losses significantly complicates the problem, since the 
seed emission itself depends on the distribution of electrons. However, Eq. (5) can be solved independently 
of the radiation field structure if SSC losses are neglected. For constant magnetic field, the solution is then 

7™a.(At) ^ , ^ ^— ^ ,^ (6) 

1 + "fmaxM/tl 

and 

N{^)^N,a-'{i-i^t/hY-\ (7) 

where At — t — z/v is the time at position z since the shock wave front passed it and v is the speed of 
the forward shock. The darkest shade of gray in Fig. 4 corresponds to At = since it is located at the 
current position of the shock wave front, where the electron energy distribution extends up to ^max- For 
times and positions within the source corresponding to At > 0, the shape of the distribution is different 
from the original if s 7^ 2 and the maximum Lorcntz factor of the distribution is less than ^max ■ The value 
of 7maa;(At) bcglns to deviate significantly from the initial value when At ~ ti/jmax, which is about 10^ s 
for B = IG and jmax = 10^. For example, for a source of size ~ 0.01 pc (light crossing time ~ 10® s) only 
the subset of electrons with Lorentz factors 7 <C 10'^ may be considered homogeneous across the source. For 
larger source dimensions the effect of inhomogeneity is even more severe. 

Light travel effects modify the appearance of the excitation structure and, therefore, the observed flux. 
It is natural to consider the internal structure in the rest frame of the shocked plasma, since it is the most 
natural frame for calculating the emission. All quantities discussed in this section are defined in the plasma 
rest frame. 

The orientation of the apparent shock fronts and the internal structure of the source for an observer in 
the plasma rest frame who is located in the y-z plane at an angle 9 to the z-axis are presented in Fig. 5. 
The equation that defines the position of the apparent shock front is 

(1 - vcosdobs/c)z = (y- R)v sin 9 obs/c + vtobs, (8) 

where tobs is the time of observation. We define tabs = as the time when the first radiation from the source 
is detected by the observer. There is no direct relation between tobs and the local emission time t because at 
any given time tobs the observer detects radiation from various locations in the source produced at different 
times. As illustrated in Fig. 5, the radiation coming from the far side of source (top of the sketch) is delayed 
with respect to that produced at the near side of the source owing to light travel time. 

The apparent shock position determines which portion of the source will be visible to the observer at 
the time of observation tobs ■ Every point in the source behind the apparent shock front contributes to the 
observed radiation at time tobs but at different local times t. A range of physical shock positions in the 
source corresponds to a single time tobs- For any given position (x, y, z) and time tobs the local emission time 
t is given by 

t = tobs + [{y - R)s[n6obs + zcos9obs]/c. (9) 

The last two terms in this equation determine the shift along the line of sight of the local point (x,y,z) 
with respect to (0, R, 0), which by definition has no delay. Once the local time of emission is found, we can 
determine At and the local distribution function of electrons, which is then used to calculate the emissivity. 
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The apparent crossing time tac is the time when the shock wave front appears to have left the source 
and is found by substituting z ^ H and y = —R into Eq. (8). This gives 

tac = 2i?sin e^bs/c + H(l/v - cos 0obs/c). (10) 

Since the speed of the shock front in the shocked plasma frame cannot exceed c/\/3, the apparent crossing 
time in the plasma rest frame can never be significantly smaller than the physical crossing time H/v despite 
relativistic speeds of colliding shocks. Higher energy electrons with fast decay times occupy a relatively thin 
region behind the apparent shock front and, therefore, the emission that they produce traces the geometry 
of the excitation region along the path of the apparent shock front. The corresponding flare peaks at half 
the apparent crossing time iac/2. On the other hand, lower energy electrons completely fill up the portion of 
the source that has been excited by the shock. Therefore, these electrons will produce emission that traces 
the entire volume of source and the associated fiare will peak at the crossing time tac- 

Synchrotron emission from the forward and reverse regions of the source can be calculated independently 
of each other at frequencies above the synchrotron self-absorption frequency. If the angle of observation 9obs 
is different from 90°, then one of the regions, reverse or forward, will block parts of the adjacent region in 
the optically thick case. 



4.2. Light Travel Time Effects for SSC Emission 

In the SSC model, a population of relativistic electrons that emit synchrotron radiation also up-scatters 
some of these synchrotron photons, thereby giving rise to inverse Compton emission. Electrons at every local 
position within the source scatter the photons that are produced locally as well as those arriving from other 
portions of the source. 

Because of this expenditure of energy, the average and maximum energy of the population of electrons 
decays with time. The intensity of synchrotron radiation produced locally therefore also decreases. However, 
even when the local electron population has decayed substantially the synchrotron emission arriving from 
other parts of the source may still contain plenty of high energy photons to drive inverse Compton emission. 

At an arbitrary point P = (xp,yp,zp) in the forward region the appearance of the shock wave front, 
which we call the internal apparent shock front, is non-planar. For time tp > zpjvj the distance from P to 
the internal apparent shock front in direction Q (see Fig. 6) is given by 

._ -/^P--P (U) 



Vf/c + cos 6 
for the forward shock front and 



Vrtp + Zp ^^2) 



Vr/c— COS 9 

for the reverse shock front in the adjacent region. Both reverse and forward regions will contribute syn- 
chrotron photons for the inverse Compton scattering. Thus, in general it is necessary to consider the 
incident emission originating from both region. 

To calculate the inverse Compton emission coefficient at point P one needs to find the synchrotron 
radiation field at this location, i.e., the intensity of the synchrotron emission as a function of direction and 
frequency. The synchrotron emission comprises light produced at different locations at different times within 
the excitation region defined by the apparent positions of the shock fronts. The time of emission from a 
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specific location (r, 0, 6) is t = tp —r/c^ whereas the shock passage time at this location is |z|/wr./, where z is 
the Cartesian coordinate of the location. Therefore, one determines the time interval At = t— \z\/vrj' between 
the emission and shock front passage for this location. The interval At determines the local distribution 
function, and hence the emitted synchrotron radiation, at this location. 



5. Calculation of Emission 

As described above, for calculating the flux in the plasma rest frame we approximate that the emission 
region has a cylindrical geometry and is filled with relativistic electrons excited by a shock wave passing 
through it at a constant speed. The density of relativistic electrons and magnetic field strength are taken 
to be constant, although the distribution of electrons changes according to the energy losses that they 
experience. The magnetic field is randomly oriented except for a possible small uniform component. A 
power-law energy spectrum of electrons is injected at the shock wave front. We use approximate expressions 
for the synchrotron emission coefficient as described in the following sections. Once the flux in the plasma 
rest frame is found, the observed quantities are calculated by using the transformation formulas discussed 
below. The following set of input parameters for each emission region is used: geometrical dimensions of 
the cylinder R (the same for both reverse and forward regions) and H] plasma number density n; magnetic 
field strength B] distribution of electrons specified by parameters 7mm: Imaxi and s; shock speed in the 
plasma rest frame v. These are complemented by the angle between the line of sight and jet axis Ooi,g] bulk 
plasma speed F' ; and redshift z, which defines the distance to the blazar. The output of the model is the 
fliix density as a function of frequency and time. 



5.1. Synchrotron Emission 

We use radiative transfer equation and the method of characteristics to calculate the intensity and 
observed flux given the synchrotron emission and absorption coefficients. Following Pacholczyk (1970), the 
synchrotron emission coefficient is 



J, (At) = c^B J iV(7, At)Fix) d-f, (13) 

where x = (71/ /7)^, Jv = ^Jv/(ciE) is the Lorentz factor of an electron (with pitch angle equal to the mean) 
that emits at critical frequency v, ci — 2.8 x 10^, and C3 — 1.87 x 10^^^ in cgs units. The integral is taken 
from 7mi„(At) to ^rnaxi^t) where At is the time since the shock front passage at the current location. The 
function F{x) = x J K^/^^z) dz, where i^s/s is a Bessel function of the second kind, can be replaced with 
its approximation (Pacholczyk 1970). After rejecting insignificant terms, we obtain 



MAt) « ^^ exp (- [7./7™a.(At)]2) ^^^^ 



czBNq 

The approximation is different from the exact formula by a factor of order unity that depends weakly on 
s, v^ and At. Assuming that the synchrotron emission by electrons provides the dominant cooling, we 
can substitute solution (6) into the above equation. To facilitate the ease of analytical considerations the 
resulting expression can be further simplified, which gives the final expression used in calculations and 

analytical considerations: 

C3_BiVo 
j^(At) =: — j^:j- exp {-vjvraax) exp {-Atjty), (15) 
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where ly'max — ci-B7,^„„j. and t^ — ti/j^^ is the characteristic decay time of electrons owing to synchrotron 
losses. A similar analysis can be carried out for the absorption coefficient, with a final result 



C3BN0 
27«meJ/2 



a„{At) = ^_ ^_ _ ^ exp {~iy/iymax) exp {-At/U), (16) 



where me is the mass of electron. The ratio of the coefficients gives the source function, which is independent 
of time and position: 

S^ = jvlajj = l^i^mev"^ . (17) 

This allows us to derive an analytical solution to the radiative transfer of the synchrotron emission. 



5.2. Radiative Transfer 

The flux density measured by an observer at distance Dq in the plasma rest frame is 



F,{tohs) = ^ Utobs,^,v)d^drj. (18) 

It is evaluated by integrating intensity observed at time tots over the projection of the source on the plane of 

the sky. The intensity along a given path defined by (^, rj) is found by solving the time-dependent radiative 

transfer equation, 

1 dl dl 

-^ + 7^^-^-^i'^ii^0+MtX), (19) 

c ot oC, 
where C is the coordinate along the line of sight. The transfer equation is integrated along a path through 
the excited region determined by the position of the apparent shock front at time tohs and the geometry of 
the source. Using the characteristics i(C) = tots + C/c, which describe a set of space-time positions through 
the source that produce the emission arriving at the observer at time tots 7 we can convert partial differential 
equation (19) into an ordinary differential equation that can be integrated with the above approximations 
for the emission/absorption coefficient. The result is 

lAtobs,^, V) = S,{1 - exp (-r,)), (20) 

with the optical depth defined by 

T,{tobs,X,v) - a.(0)ACT,(x)cxp(-At/t,). (21) 

Here, r^(x) = (expx - l)/x and 

X - {cosOobs - v/c)AC/{vU). (22) 

The length AC = (max — Cmin is the geometric thickness of the source along the given path through it. By 
combining Eqs. 18-22 with the expression for At, 

At = tabs +?7 sin 6'ofcs/w + (1/c- COs9obs/v) (rmn, (23) 

one can obtain an analytical expression for the evolution of the observed synchrotron spectrum if the geometry 
of the source is sufficiently simple. For a cylindrical source, the expression for AC is greatly simplified if the 
direction toward the observer coincides with the axis of the jet. In this case the solution is 

FAtobs) = ^[1 - exp(-r,)]5„ (24) 
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with 

T^ = ai.(0)i?exp(-tofcs/i,,) — , (25) 



where tac ~ (c/v — 1)H / c for the forward region. 



5.3. Self-Compton Emission 

Similar to the synchrotron case, the inverse Compton flux density is obtained by integrating intensity 
over the extent of the source on the plane of the sky. Since we assume that the inverse Compton emission 
is optically thin at all frequencies of interest, the intensity is the integral of the inverse Compton emission 
coefficient along a given path toward the observer, which gives 

F^itobs) = ^J J Jj?(tA,V,Od^dvdC, (26) 

where t{() = tabs + C/c takes into account the light travel effects. 
The inverse Compton emission coefficient 

j^ = fdn, f djN{j) I dv, ^i,^{fi,)^b{iy,,i^,n,,nf,j) (27) 

depends on the local distribution of electrons N{"/), which is assumed to be isotropic, and the incident 
synchrotron radiation /jy. (O^) that provides the seed photons at different frequencies Vi. The redistribution 
function tjj determines the fraction of incident radiation at initial frequency Vi in initial direction f2i that is 
scattered by electrons with Lorentz factor 7 into final direction Uf with final frequency ly. The final direction 
rif is constant and determined by the angle Oobs- The limits of the integrals, as well as the expression for 
the redistribution function, are discussed in Appendix B. The redistribution function that we use applies 
at all electron energies and photon frequencies allowed by energy and momentum conservation. It does not 
make any simplifying assumptions about the angular distribution of the incident radiation, and is accurate 
in both Thomson and Klein-Nishina regimes. li jhi^i <C nieC^ (Thomson limit) and Vi <C i^f, expression (27) 
for the Compton emission coefficient simplifies to that derived by Reynolds (1982). 

Similar to the calculation of the observed emission, the internal synchrotron seed photon intensity is 
calculated by solving the time-dependent radiative transfer equation for each direction fli. The characteristic 
equation is t(r) = t — r/c^ where r is the distance from the point at which the intensity is evaluated in the 
direction opposite to fJi, i.e., the direction from which the radiation comes. If the same approximations for 
the absorption and emission coefficients are employed, the resulting intensity is 

I,Xt,x,y,zA) = [l-exp(-T,J]5,., (28) 

where 

T,, =a,,(0)L(a)T..(x)exp(-At/t,J (29) 

and, for the incident emission from the forward region, x = (^//c — ^iz) L f (fli) / (v fti,.) and At = t — z/vf. 
Here, Lf{fli) is the distance to the edge of the excited material in the forward region and fliz is the 
projection of vector fli onto the z-axis. For the incident emission originating in the reverse region, x = 
{vr/c + fliz)Lr{i~ii) / {vrt^.) and At ^ t — z/(cCliz). 

In calculating the inverse Compton flux density we use the analytical solution for the incident syn- 
chrotron intensity and perform all calculations on an adaptive grid that takes into account the electron 
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excitation structure and the internal synchrotron emission field, which can be highly inhoniogeneous and 
anisotropic. 



5.4. Observed quantities 

The model calculations are performed in the plasma rest frame where a fiducial observer is located at 
distance Dq along a line that subtends an angle 9obs to the z-axis. The flux density thus calculated must 
be converted to that received in the observer rest frame (denoted by star) by taking into account all the 
relativistic, light travel, and cosmological effects. 

The relation between the angle of observation in the plasma rest frame 9obs and in the observer frame 
^obs is 

vL/c + COS^ohs 

cos^:,,--^^^ ^. (30) 

1 + v'p cos Oobs/C 

For 6obs = 90°, the viewing angle in the observer rest frame can be found from cos0*j,^ = v' /c, which gives 
^obs ~ l/^p! this is the angle that maximizes apparent superluminal motion for a fixed speed v' . 

The transformation formulas for frequencies, time, and flux are simplified with the use of the Doppler 
factor, defined as 

^^T'^ii-v'^coseijcy ^^^) 

Thus defined, the Doppler factor, combined with cosmological time dilation, gives the following relations: 

z.* = lyS/il + z), Cb. = (1 + z)tobs/6, (32) 

and 

F: = {Do/D*fd''F,, (33) 

where D* is the luminosity distance to the blazar measured by the observer. Accordingly, the apparent 
crossing time t*^, = (1 + z)tac/6 in the frame of the observer. 



6. Model Light Curves and Spectra 

In this section we present the result of our simulations. We first consider the light curves and spectra 
from a single region of emission. Contributions from the adjacent region can be neglected if it is geometrically 
thinner and/or contains lower density plasma and weak magnetic field. We compare synchrotron and inverse 
Compton flares with respect to possible time delays between high and low energy emission. We finish this 
section by presenting results of a double region simulation and our modeling of variability from the quasar 
3C 273. All the quantities discussed here, except for the flare in 3C 273 will be given in the rest frame of 
the emitting plasma. Corresponding quantities in the observer's frame can be easily calculated if the speed 
of the plasma, orientation of the jet, and source redshift are known. 



6.1. Synchrotron: Single emission region 

The evolution of the synchrotron spectrum during a flare is shown in Figs. 7 and 8 for two representative 
viewing angles. The shape of the synchrotron spectrum and its evolution can be roughly described by 
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examining three critical frequencies: 

1. The synchrotron self-absorption frequency i>t, at which the source becomes optically thick, is found by 
solving the equation ay(0)i = 1, where L is the average geometrical thickness of the source along the 
line of sight. The solution is 



fcsNoL 



s + 4 



s + 2 



"'^[^r^ i^^B) — . (34) 

During the flare the thickness L is expected to grow as a result of the shock wave propagation through 
the excitation region. Even when the viewing angle 9obs = 90°, light travel effects cause the apparent 
shock wave to propagate along the line of sight into, as well as across, the excitation region. The 
growing L will result in the increase of the self-absorption frequency vt until the whole excitation 
region has been traversed by the shock. 

2. The expression for the break frequency Vb at which the spectrum steepens owing to radiative energy 
losses depends strongly on the viewing angle. Analytical solution (25) of the transfer equation for 
Qobs = 0° allows one to find the break frequency by solving the equation t,^ = jmTi{tobsitac}, which for 
tabs < tac gives 

z/b=cii?/(iLB). (35) 

Substituting tabs = tac and the input parameters used in simulations (see Fig. 7) gives Ub = 6.3 x 10^^ Hz. 
At the beginning of the flare the spectral index, defined by F^{tobs) oc z^", is —0.5 at all frequencies 
between the self- absorption and the drop-off frequencies but it decreases quickly to —1 as the break 
frequency sweeps across the spectrum until it reaches its minimum value of 6.3 x 10^^ Hz at tobs = tac- 
The spectral index steepens by 1/2 above the break frequency due to the fact that the depth of the 
actual emitting volume ai v > Vb is determined by the decay time ti, ex v^^'"^ . Similar reasoning may 
be used for different viewing angles, in which case Eq. (10) should be used to find the appropriate 
crossing time tac- We stress that the break in our simulated spectrum is only due to the frequency 
stratification of the emitting electrons. 

3. The exponential drop-off frequency Vd remains equal to the critical frequency of the maximum-energy 
electrons, Vmax = (^iB^maxi as long as the apparent shock wave is in the excitation region. However, 
after the shock exits the excitation region the decrease in the drop-off frequency mimics the decrease 
in the maximum electron energy (cf. Eq. (6) with At = i — tac)- 

The assumption that the excitation region is a cylinder oriented along the axis of the jet affects the 
behavior of the observed spectrum as a function of the viewing angle. At viewing angle 9obs = 0° the 
apparent surface of the shock propagates along the axis of the cylinder and, therefore, the circular geometry 
of the source has no effect on the observed spectrum. The evolution of the spectrum simply reflects the 
increase in depth through the source as the shock wave propagates through it. On the other hand, for 
viewing angle Oobs — 90° the apparent surface of the shock front is tilted away from the plane normal to 
the cylinder axis. Thus, the apparent surface traces a roughly circular shape for such a viewing angle; this 
causes significant differences in the evolution of the observed spectrum relative to the Oobs = 0° case. In 
particular, emission at the highest frequencies peaks roughly at half the apparent crossing time, since the 
intersection between the excitation region and the shock front is largest at about that time. The flux density 
will decrease significantly at these frequencies by the time the shock front exits the source. At the crossing 
time tobs — tac, the spectral index assumes a value of about —1.75 at the highest frequencies owing to a 
unique combination of the cylindrical geometry and the orientation of the excitation structure in the source. 
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This stage is very brief and requires exact cylindrical geometry; hence, it is unlikely to be observed but can 
be used to validate numerical calculations. 

The spectral behavior is easier to interpret if one compares it with a set of complementary light curves 
at a few representative frequencies. At zero viewing angle (Fig. 9), all light curves peak at the crossing time 
tac ~ 8 X 10^ s for the parameters used (cf. Fig. 7). However, only those light curves for which the decay 
time is close to the crossing time, t^ ^ tac, will be roughly symmetric and have a sharp peak; i.e., one should 
expect symmetric flares only when the frequency of observation is close to the break frequency at io^s = tac 
for zero viewing angle. At higher frequencies v ^ vi,, one will observe a sharp rise followed by a flat top and 
then a quick decay. At lower frequencies i^ <C i^b, an increase in flux density that occurs on time scale tac will 
be followed by a relatively long quasi-exponential decay on time scale t^. In the latter case, one will observe 
the emission rise to a relatively stable high value of flux with insignificant decay if the period of observation 
is comparable to tac- Although there is no time delay in the peak emission in the case of zero viewing angle, 
a cross correlation analysis should detect a positive time delay between high and low frequency synchrotron 
flares. This is because the flares at lower frequencies have larger decay times and do not respond to changes 
in the source as quickly as the flares at higher frequencies; i.e., lower frequency flares both rise and decay 
more slowly. 

For viewing angle Oobs — 90° (Fig. 10) both the shapes and the positions of the peak emission are different 
from the zero viewing angle case. Due to the circular shape of the excitation region and the orientation of 
the apparent shock front, the flares appear symmetric even at frequencies higher than the break frequency 
if the source is geometrically thin along the axis of the jet {H <C R). The peaks of the flares occur between 
half the crossing time iac/2 and the crossing time tac, with higher frequency flare peaks tending toward ioc/2 
and those at lower frequencies approaching tac- The time delay of the peak emission is never expected to 
exceed half the crossing time or roughly half the duration of the flare. While comparing Figs. 9 and 10 one 
should bear in mind that the crossing time tac = 2 x 10^ s for viewing angle 9obs — 90° is longer by 2.5 than 
that for zero viewing angle. 



6.2. Inverse Compton: Single emission region 

The inverse Compton (IC) emission coefficient in the SSC model results from the convolution of the 
electron energy spectrum and the synchrotron emission spectrum. The incident synchrotron intensity may 
depend strongly on the direction of the incident emission fi^, with each direction having its own set of char- 
acteristic frequencies (f^, i^t,, i/^) for a specific time and location in the source. However, the self-absorption 
frequency is not expected to depend as strongly on direction or time and, therefore, its average value may 
be used in a qualitative analysis. The turn-over frequency of the IC emission can therefore be estimated as 
Vet — '^lmin{^t)i whcrc {vtj \s the average self-absorption frequency of the internal synchrotron radiation 
field. The synchrotron drop-off frequency Vd for the directions terminated by the internal apparent shock 
front is constant and equal to the critical frequency of highest energy electrons, Vmax = ciB^f^ax- Thus, the 
IC emission is expected to drop off exponentially at frequencies greater than Vc^max = '^Imax^max as long as 
the internal shock front is inside the excitation region. Once the internal shock front exits, the IC drop-off 
frequency Vc^d begins to decrease. The break frequency of the incident synchrotron emission v^ is a strong 
function of direction and time and thus does not translate into a well-defined break frequency of the IC 
spectrum. Rather, we expect the IC spectrum to steepen gradually by Aa « 1/2 toward higher frequencies. 

The evolution of the IC spectrum for two different viewing angles is shown in Figs. 11 and 12. The 
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overall behavior is similar to that of the synchrotron emission. However, instead of the well-defined power- 
law seen in the synchrotron spectrum the IC spectrum gradually steepens over a wide range of frequencies. 
Hence, it cannot be approximated as a single power-law over more than two decades in frequency. It is also 
clear that the IC emission peaks after the crossing time at lower frequencies (see dashed line in Fig. 11). 
This pattern is more obvious if one examines a set of IC light curves at different frequencies, as shown in 
Figs. 13 and 14. For viewing angle 6*06^ = 90° (Fig. 14), the behavior of the IC light curves is very similar 
to that of the synchrotron flares. The emission at higher frequencies tends to peak at time ^ tac/2. The 
peak shifts toward the crossing time tac at lower frequencies. The light curves appear symmetric (similar 
rise and decay time) over a broad range of frequencies. This suggests that the internal light travel effects 
are insignificant for this viewing angle and that the shapes of the flares are determined by the geometry 
of the excitation region along the path of the apparent shock front. However, in the case of zero viewing 
angle (Fig. 13) the results of our simulations demonstrate that the IC flares at sufficiently low frequencies 
peak after the crossing time, in contrast to the synchrotron flares at this viewing angle, which all peak at 
the crossing time. Due to internal light travel delays, the IC flares can still increase in brightness after the 
apparent shock front has left the excitation region and the high energy electrons are no longer supplied. 

The decay time of the IC emission is difficult to define uniquely because it depends on the frequency- 
integrated incident synchrotron radiation, whose characteristics may depend strongly on direction, time, and 
location in the source. Even when the position of the shock front as seen by a distant observer appears to 
have left the source, the internal shock front as viewed from a given location may still lie well inside the 
excitation region for a large range of directions, thereby supplying synchrotron seed photons at the highest 
frequencies. An upper estimate of the IC decay time can be obtained by assuming that the maximum 
incident synchrotron frequency Vmax does not depend on time. Then for a given IC frequency Vc the decay 
time is the time needed for the drop-off frequency of the inverse Compton emission Vc^d = ^li^tYvmax to 
decrease to Vc- This gives the estimate 

i^^^t /l^. (36) 

V ^c 

The actual decay time may be shorter because of the decrease in the drop-off frequency of the incident 
synchrotron radiation. 



6.3. Time delays 

Comparison of the synchrotron and inverse Compton light curves offers a possible explanation of the 
existence of negative time delays between radio and X-ray flares. In rough terms, the time delays depend on 

(1) the apparent crossing time tac^ which is determined by the source dimensions and the viewing angle, and 

(2) the decay times for synchrotron and IC emission. If both decay times are much smaller than the crossing 
time, then no significant shifts between these light curves will be present. Decreasing the frequency of 
observation of the synchrotron emission will lead to an increase in the decay time. If it becomes comparable 
to or longer than the crossing time, the synchrotron fiare will be delayed with respect to the IC fiare. If, on 
the other hand, one compares synchrotron emission at sufficiently high frequencies (infrared or optical) with 
IC emission at relatively low frequencies (soft X-rays) so that the corresponding IC decay time is longer than 
the crossing time, then the IC flare should peak after the synchrotron flare. 
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6.4. Effects of the dual emission regions 

As a result of the collision between shocks both reverse and forward regions will be created. Thus, the 
observer may see a combination of emission from these two regions. For the most energetic collisions both 
shock speeds should be close to c/3 in the rest frame of shocked plasma. Within the assumed cylindrical 
geometry, one naturally expects the cylinder radius to be the same for both regions, R = Rr = Rf, but the 
heights Hr and Hf may be different. 

Due to intrinsic geometrical symmetries in the problem, one can conclude that for viewing angle 6obs = 
90° the synchrotron emission observed from both regions will be identical if all the input parameters of both 
regions are the same, i.e., if Hr = Hf, rir = nj, B,. = Bf, etc. For all other viewing angles, if the apparent 
crossing times for both regions are equal, tac.r = tacj, which imposes a constraint on the values of Hr and 
Hf, and if the other input parameters of both regions are identical (i.e., rir — Uf, Br ~ Bf, etc.), then the 
optically thin emission from both regions will have identical shapes, albeit the flux level should be different 
because of the difference between Hr and Hf. By adjusting A^o = («— l)?^7^^/[l— ilmin/lmaxy~^] one can 
change the peak flux levels without altering the shape of the light curve. On the other hand, even moderate 
discrepancies in s, jrnax, or B may lead to dramatic differences in the profiles of the reverse and forward 
flares. 

If observed flares in blazars stem from collisions between shocks, then all the flares must consist of two 
components (reverse and forward). In most cases one of the components will overwhelm its counterpart, 
resulting in an apparently single profile. Furthermore, since blazar flares occur frequently, it is challenging 
to separate doubly peaked flares resulting from a single collision and two separate events occurring at about 
the same time. 

To illustrate the possible morphology of a doubly peaked flare, we have conducted a simulation with the 
forward region being denser but thinner than the reverse region (Fig. 15). Synchrotron seed photons origi- 
nating from the adjacent region were taken into account in calculation of the IC light curves. Examination 
of the results indicates that the seed photons in the reverse region are dominated by synchrotron emission 
from the adjacent forward region because of the higher density and magnetic field in the latter. The light 
travel delays for the seed photons from the adjacent region are also longer on average. The choice of the 
input parameters results in a quick spike of synchrotron emission (forward region) followed by a weaker but 
longer outburst (reverse flare). The corresponding IC variability has a more complex profile because the 
IC emission from the reverse region is strongly affected by seed photons from the adjacent forward region. 
Notice that both flares begin at the same time, in accordance with our model. The shock collision model 
cannot explain a doubly peaked flare consisting of a broad weak outburst followed by a spike of emission at 
the end because that would require different onset times for the reverse and forward flares. 



6.5. Application: 3C 273 

In March 1999 the quasar 3C 273 produced a short outburst that was observed simultaneously in the 
infrared K-band and the X-ray band. An examination of the multifrequency data (Fig. 16) shows that the 
X-ray flare was delayed with respect to its K-band counterpart by about one day. The activity in K-band 
spanned about two days, whereas the X-ray flare continued for about four days. The spectrum of the X-ray 
emission was observed to steepen during the flare. 

This long time delay can be explained by applying the shock collision scenario, according to which the 
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collision leads to the formation of two adjacent emission zones with similar properties (reverse and forward 
regions). Two conditions are required to produce a successful fit: (1) the X-ray flare is dominated by IC 
emission from a larger region that has a crossing time of about four days, and (2) the infrared flare is in 
turn dominated by synchrotron emission from a smaller region with apparent crossing time of about 2 days. 
In general, these conditions are difficult to satisfy simultaneously since brighter synchrotron emission will 
correspond to brighter IC emission. However, if the Lorentz factor of the maximum energy electrons is 
limited to jmax ^ 10'^, the magnetic field strength i? ~ 1 G, and the bulk Lorentz factor of the emitting 
plasma F ^ 10, then even a small discrepancy in B and/or jmax between the forward and reverse regions 
may significantly affect the K-band synchrotron emission. 

For simplicity and in order to minimize the number of freely adjustable parameters, we adopt a bulk 
Lorentz factor F = 10 and angle between the jet axis and line of sight 9obs = 90° in the plasma rest frame. 
This viewing angle corresponds to 0*{,^ = 5.7° in the observer's frame. We also assume that the speeds of 
both the forward and reverse shocks equal to 0.34c in the rest frame of the emitting plasma. The observed 
duration of the flares combined with the direction of motion and the bulk speed of the plasma sets the size 
of the emission region at about 0.5 x 10^ sec, whereas the roughly symmetric shape of the flares sets the 
magnetic field at about 1 G. Adjusting these quantities and comparing the results of the simulation with 
the data culminate in the fit presented in Fig. 16. 



7. Discussion 

7.1. The Nature of Radio Core 

The model that we have developed can be seen as part of a more general picture of the jet in which a 
stationary system of oblique shock waves terminated by a Mach disk corresponds to the stationary feature 
called the radio core, which is seen in essentially all jets. The model therefore relates the gas dynamical 
structure of the jet to its variability properties. Provided that the shock collision scenario is correct, high 
speeds of superluminal components originating from the radio core (bulk Lorentz factors of 10 and higher) 
would require relativistic shocks upstream of the stationary shock complex with bulk Lorentz factors up to 
200 or even higher, since the collision will tend to substantially decelerate the moving shock. This requirement 
of high upstream speed may be relaxed if interaction with the oblique shocks alone can energize the moving 
shock sufficiently enough to explain the superluminal features, since oblique shocks are not expected to slow 
down the moving shock appreciably and, therefore, extremely high speeds prior to collision are not necessary. 
We have only conducted a ID study study of the collision, but the system of oblique shocks and Mach disk 
is inherently a 3D structure. Therefore, more realistic simulations are required to ascertain the plausibility 
of our scenario and its gas dynamical consequences for relativistic jets. 

The results of the modeling reported in this paper do not depend on the details of how the excited 
plasma is produced. Most calculations are performed in the rest frame of the emitting plasma, which gives 
them considerable generality since one can transform to a frame in which the plasma is moving. Only the 
fit to the 3C 273 flare has been converted into the observer frame to facilitate easy comparison with the 
observed light curves. On the other hand, we have applied the model to a scenario in which the flares are 
caused by a collision between a moving shock generated at the base of the jet and a stationary Mach disk 
located downstream in the jet at a distance of at least ^ 1 pc. This collision scenario relies on the existence 
of a stationary shock wave complex, and hence is different from the "internal shock" scenario investigated 
by Spada et al. (2001), who considered collisions between moving shells produced at the base of the jet. The 
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moving shock in our scenario and the shcUs in the model of Spada et aL (2001) could both be produced at 
the base of the jet by similar mechanisms, perhaps involving an instability in the accretion disk. A shell of 
relativistic plasma is formed behind the moving shock as it propagates downstream in the jet. The collision 
of the moving shell with the Mach disk will occur at roughly the same location in the jet, i.e., the location 
of the Mach disk, regardless of the frequency of ejections at the base of the jet. In contrast, a collision of 
two moving shells will occur at different locations in the jet, since the time required for the faster shell to 
catch up with the slower shell depends on the relative speeds and initial separation of the shells. However, 
it can be easily shown that the ejection of two shells separated by ^ 5 days will result in the collision at the 
distance of ^ 1 pc if the Lorentz factor of the slower shell is ~ 10. Because of this, it might be difficult to 
distinguish observationally between these two scenarios. Furthermore, our formalism can be applied to the 
Spada et al. (2001) scenario with a suitable Lorentz transformation. 

The bright appearance of the superluminal knots seen in VLBI images is at least partly due to the 
relativistic speed of the plasma producing the emission. In the majority of blazars the radio core usually 
appears brighter than the superluminal features. However, a stationary shock wave normal to the flow of 
plasma in the jet slows the flow down to a mildly relativistic speed. Therefore, minimal relativistic boosting 
is expected in the case of quiescent emission from the plasma on the downstream side of the Mach disk. The 
radio core is likely to be dominated by relativistic flow passing through the system of oblique shocks or the 
reaccelerated flow farther downstream of the Mach disk. 

The outcome of the collision depends on the relative speeds of the moving shock and the quiescent flow 
between the shock and the Mach disk prior to collision. If the speed of the moving shock prior to collision 
is much higher than the speed of the quiescent flow, then the moving shock will destroy the stationary 
shock complex, resulting in an excited plasma propagating downstream at relativistic speed. This scenario 
requires that the oblique shocks and the Mach disk reestablish at approximately the same locations in the 
jet after the collision. The time for this system to reappear in the jet is quite short (see Gomez et al. 1997). 
The collision will manifest itself as an increase in brightness of the radio core followed by the appearance 
of a superluminal knot associated with the shocked plasma that drags the core downstream. This must be 
followed by the reappearance of a bright feature at the original position of the radio core, signifying the 
reestablishment of the Mach disk. It might be possible to detect the temporary downstream translation of 
the apparent core in high resolution, multi-epoch VLBI observations with phase-referencing to register the 
coordinates of the images. Alternatively, in the case of a weak moving shock, the stationary shock complex 
will be weakly affected by the collision. It will retain its position in the jet with probably only a slight shift 
downstream. In this case, a superluminal knot associated with the flare can be identified with the portion of 
the moving shock that only interacts with oblique shocks between the Mach disk and the outer boundaries 
of the jet. 

The existence and properties of the Mach disk established in 2D simulations and analytical studies must 
be confirmed with high resolution 3D simulations. We have assumed in our studies that the oblique shocks 
play an insignificant role in the observed variability. However, existing 3D simulations indicate that strong 
oblique shocks can appear and that the Mach disk may be small or even non-existent if fluctuations of the 
jet direction are introduced (Aloy et al. 2003). This and other issues such as the time scale on which the 
stationary structure can be reestablished should be addressed by future jet simulations. 
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7.2. Emission calculation 

Our simulations rely on approximations to the synchrotron emission and absorption coeffieients that 
offer a number of advantages. They allow a more immediate yet comprehensive description of the excitation 
properties of the medium energized by a propagating shock wave. We have used these approximations to 
derive analytical solutions to the time-dependent radiative transfer equations that describe the propagation 
of synchrotron emission through an evolving relativistic medium with all light travel time effects taken into 
account. This significantly speeds up the fC emission calculation. For the cases involving simple geometry, 
such as a cylindrical source at zero viewing angle, the equations can be integrated analytically. Despite these 
advantages, the approximations can give erroneous results at the later stages of the flare when the maximum 
electron energy ^maxi^t) ~^ Imin- To simplify the analytical solutions, we have substituted exp(— [Ai/ij^]^) 
with exp(— Ai/ij^) to make the IC emission calculations with light travel delays for the seed photons possible. 
This procedure can underestimate the rate of decay of the observed synchrotron emission, especially at low 
frequencies. A more precise formulation can be worked out based on direct integration of exp(— [Ai/i^]^) 
and expression of the results in term of the Error function. 

We have assumed that the electron energy losses are dominated by synchrotron emission. Although this 
may not be the case in some sources, incorporation of inverse Compton losses in a calculation that relies 
on a dynamical evaluation of the time-delayed synchrotron emission poses considerable difficulties. Indeed, 
if inverse Compton losses dominate then the electron distribution function at a given location and time in 
the source will depend on the synchrotron radiation field at this location. This photon field in turn depends 
on the structure of the electron energy distribution throughout the source at retarded times, as needed to 
account for the light travel delays. Thus, the problem of determining the distribution of electrons not only 
becomes non-linear but also non-local, since the distribution of electrons at a given position will depend 
on the electron excitation structure in the whole emission region at different retarded times. This task lies 
beyond the scope of the current work. 

Source expansion, adiabatic energy losses, and gradients in magnetic field and electron density have 
been excluded from the current work as well. Consideration of these effects must be postponed until a more 
realistic 3D picture of the collision becomes available. The IC emission calculation relies on a number of 
analytical estimators that apply to a non-expanding source with uniform magnetic field and density, a shock 
wave moving at constant speed, and evolving electron energy distribution dominated by synchrotron energy 
losses. Many of these estimators must be reconsidered if the source is allowed to expand substantially during 
the flare, or if gradients in magnetic field or electron density are introduced. 

In this paper we have concentrated on the SSC mechanism for production of X-rays. External Radiation 
Compton (ERC) emission, which is not affected by the light travel time effects in the same manner as in 
the SSC model, will be considered in a separate work. Since external emission in most cases is considered 
constant, it is expected to fill the source uniformly. Despite the absence of internal light travel delays 
in the ERC model, one cannot reject this model based solely on the presence of the relative time delays 
between multifrequency fiares: the internal excitation structure due to electron energy losses can produce 
both positive and negative time delays between synchrotron and IC variations for certain viewing angles, 
most notably Oobs ~ 90° in the plasma rest frame. However, for viewing angles close to 6obs =■ 0°, there is a 
fundamental distinction between SSC and ERC models. The ERC emission is expected to peak at crossing 
time tac similarly to the synchrotron flares at this viewing angle. The SSC emission, on the other hand, can 
continue to brighten even after the crossing time because of the light travel delays of the internal synchrotron 
emission. Mirror Compton emission, even if strong enough to compete with SSC emission, is expected to be 
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delayed by at least several times tac and is, therefore, completely disconnected from synchrotron and SSC 
variability. 

A preliminary investigation indicates that the 3C 273 flare discussed above cannot be explained by the 
ERC model with emission from the dusty torus as a dominant contribution. The associated decay times tend 
to be higher than in the SSC model because the maximum frequency of available seed photons is higher and 
does not decline with time, resulting in flares with very long decay time, in contrast with observations. The 
inclusion of adiabatic losses will decrease the time scale of variability of the ERC flare. More importantly, 
the change in Doppler boosting of the external radiation as the emission region propagates farther along the 
jet may quench the ERC flare before any substantial decay of the energy of electron producing the emission. 
These two effects must be significant for the flares of ^ 1 month or longer in duration. However, their effect 
on the fast flares such as the one in 3C 273 that we have considered may be neglected. 



7.3. Comparison w^ith Observations 

Light travel delays strongly influence the multifrequcncy timing of variable emission produced in rcla- 
tivistic jets of blazars. Despite considerable complications that the light travel delays present in modeling 
emission variability, they offer a number of diagnostics of the emission regions that can be used in interpreting 
observational data. In this section, we summarize some of these diagnostics. 

Variability of synchrotron radiation at high frequencies originates from a thin emission zone around the 
shock front propagating through the excitation region. Thus, this emission traces the geometrical shape of 
the excitation region along the path of the apparent shock front, whose orientation as a function of viewing 
angle is given by Eq. (9). For instance, for a viewing angle of 90° in the plasma rest frame, which may 
correspond to only a few degrees in the observer's frame if the jet is highly relativistic, the apparent shock 
front is tilted with respect to the orientation of the physical shock front (see Figs. 4 and 5). Because of this 
tilt, the emission does not simply trace the geometry along the jet, but instead depends on the shape of the 
excitation region both along and across the jet. If the excitation region is thin [H <C i?), the shape of the 
observed flare will reflect its structure across the jet. The geometry along the axis of the jet can only be 
probed in isolation if the viewing angle is very small in the plasma rest frame, which requires almost perfect 
alignment of the jet axis and the line of sight in the observer's frame. 

Negative time delays between synchrotron and IC flares, with variability at lower frequencies leading 
those at higher frequencies, can be used to constrain the emission models. As our simulations indicate, for 
relatively small viewing angles the internal light travel delays may exceed the apparent crossing time. In this 
case, the IC flux at sufficiently low frequencies will continue to increase for some period after the crossing 
time. On the other hand, the ERC model is not expected to exhibit such behavior. To estimate the range of 
viewing angles consistent with such behavior we plot in Fig. 17 the apparent crossing time tac as a function 
of viewing angle and an estimate of the internal light travel delays based on the maximum time delay for 
the central point in the source (0, 0, H/2): 



Unt - H/v + yfW+{H/W/ c. (37) 

For a thin source {H = 0.2R) the viewing angle in the observer frame must be significantly smaller than 
^ 35°/(2r' ), where F' is the bulk speed of the emitting plasma. 

For viewing angles 9oi,s ^ 90°, the internal light travel delays are relatively small compared to the 
crossing time tac- For such viewing angles, frequency stratification and geometry of the source will determine 
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relative time delays between flares at different frequencies. In general, the ERC flares may exhibit similar 
behavior to the SSC flares since the electron excitation structure and source geometry are the same. However, 
the ERC model can still be rejected as a viable X-ray emission model for different reasons. For example, if 
the external emission is provided by dust particles in molecular torus at temperature Td = 1000 K, the seed 
emission will be dominated by infrared photons. In this case, low energy relativistic electrons (7 ^ 15) are 
sufficient to produce X-ray emission in the RXTE band. These electrons decay very slowly and, therefore, the 
ERC flares observed with RXTE are expected to have extremely long decay times, rendering them incapable 
of fitting symmetric flares with roughly similar rise and decay times as seen in 3C 273 (see Fig. 16). 

The internal light travel delays are emphasized even further when both reverse and forward zones 
are taken into account, since some seed photons should arrive from the adjacent region, which will require 
additional time. Because of the seed emission originating in the adjacent zone, IC flares may possess complex 
profiles that cannot be described as a sum of two independent flares in the same way as can be done with the 
synchrotron flares. Moreover, as our modeling of the 3C 273 flares indicates, longer observed time delays can 
be explained with SSC emission from dual regions than is allowed by single region modeling. In particular, 
it is possible to explain SSC flares that peak not only after the maximum of the synchrotron emission, but 
also when the synchrotron emission has virtually returned to the background value that existed prior to the 
flare. Further investigation of SSC flare profiles from dual emission regions and their dependence on viewing 
angle will be considered in a later paper. 



8. Summary 

We have developed a model of the emission regions in blazar jets responsible for the fast outbursts 
of flux that characterize this class of active galactic nuclei. The collision of a moving shock wave and a 
stationary shock complex in the relativistic jet results in a pair of propagating shocks that drive acceleration 
of particles and strengthen the magnetic field, thus producing enhanced synchrotron and IC emission. The 
model incorporates time-dependent radiative transfer through a rapidly evolving relativistic medium to find 
the time-delayed synchrotron radiation field that provides seed photons for IC scattering. The ability of the 
model to trace the evolution of the excitation structure in the source provides a powerful tool for probing the 
inner geometry of the jet structures responsible for the sudden outbursts of activity usually associated with 
the ejection of superluminal features from the radio core. The model is especially well suited for modeling 
multifrequency variability involving relative time delays between light curves in different spectral bands. 

The investigation of internal light travel delays that apply to the seed emission in the SSC model strongly 
suggests that the relative time delays between synchrotron and inverse Compton flares can be dominated by 
this effect. The double structure of the excitation region resulting from a shock collision can explain extreme 
reverse time delays when the peak of the inverse Compton emission occurs after the synchrotron flare, which 
starts at the same time, has returned to the background level. 

This research was supported through NASA grant NAG5-11811 and National Science Foundation grant 
AST-0098579. 
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A. The Apparent Position of the Shock Front 

Expression (8) for the apparent position of the excitation front may be derived as follows. We define 
coordinate system (^, ry, C) and consider emission from an arbitrary point P = {x, y, z) within the boundaries 
of the source (see Fig. 18). The coordinate system (^,?7, C) is obtained from the (x,y,z) coordinates by 
shifting the origin to (0, R, 0) and rotating the shifted z-axis by angle 9obs in the x = plane such that the 
(^ axis points toward a distant observer in the plasma rest frame. 

The shock front reaches point P at time ti — z/v and excites the electrons there. A time interval 
^2 = [{R — y) sin 9obs — z cos 9obs]/c later, the radiation produced at P reaches the C = plane that intersects 
point (0, R, 0) and is normal to the line of sight. A further time interval ts = Dq/c is required for this 
radiation to reach the observer at distance Do from the plane C = 0. Thus, the observer will detect radiation 
from point P at time 

tots =ti+t2+ h - Dq/c = z/v + [{R - y) sin Oobs - z cos dobs\/c- (Al) 

The time tobs is defined such that the observer detects emission from point (0, i?, 0) at time tobs — 0. 
Rearranging terms in this equation yields 

z(l - i; cos 6*063/2) = (y- R)v'&m6obs/c + vtobs, (A2) 

which defines the apparent shock front position at time tobs ■ 



B. The Redistribution Function and Limits of Integration 

In general the limits of integration in Eq. (27) are determined by conditions under which it is energetically 
possible for an electron with a Lorentz factor 7 to scatter a photon of frequency Vi moving in direction il^ 
into new direction $1/ with new frequency Vf. For an isotropic distribution of electrons, the condition is 

(72 - l)[(e, - e/)2 + 2e,e/(l - ^ls)] > [-f{e, - e/) - e,e/(l - fi,)f , (Bl) 

where Cij- = hvij / {rrieC^) and /i^ = il^ • Vlf. 

We will consider the limits of the integration in Eq. (27) in the Thomson limit when 76^ ^ 1. This results 
in pure Thomson scattering in the rest frame of the scattering electron so that the photon frequency does not 
change in this frame of reference. The relativistic transformations from the rest frame of the plasma to that 
of the scattering electron and back to the plasma rest frame are the only factors responsible for the photon 
energy change in the Thomson limit. The highest synchrotron photon frequency available for scattering in 
the SSC model is Vi^max ^ ciBj"^^^. For 7^02: — iC and B — IG one finds i^i^max ^ 4.2 x 10^'* Hz and 
lmax&i,max ~ 0.02 ^ 1. Thus, for typical values of the magnetic field and maximum electron energy found 
in blazars, inverse Compton scattering occurs in the Thomson limit. 

Applying formula (Bl) to the particular integration order in Eq. (27) provides the following limits. 
Both outer integrals are taken over the complete range of possible values, i.e., f^^ is integrated over the entire 
sphere and 7 is integrated from 7mi„(Ai) to ^max{^t)- In practice, the lower limit on 7 may depend on the 
maximum available synchrotron frequency and the final frequency of scattering. Given the final frequency 
Vf, scattering angle /ig, and scattering electron Lorentz factor 7, the range of energetically possible values 
of the initial frequency vi is limited to 

Vf/A <Vi< Avf, (B2) 

where A = [1 + (7^ - 1)(1 - ^,)] + ^[1 + (7" - 1)(1 - l^sW - 1- 
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For /is = 1, i.e., no change in direetion of the photon, the only allowed value is Vi — Vf] no scattering can 
effectively occur in this case. On the other hand, for /is = — f the maximum possible range of the final 
frequencies is achieved, namely z^//(47^) < Vi < Aj'^Vf for 7 3> f. In this case, the largest gain of the 
photon's energy Cf — Aj'^Ci occurs in a head-on collision with an electron, while the photons that approach 
electrons from behind lose most of their energy: Cf =■ 6^/(47^). In practice, the upper limit on Vi will be 
determined by Vmax- 

Applying relativistic transformations, combined with the Klein-Nishina cross-section, to the problem of 
scattering results in a rather complicated integral expression that includes the incident radiation spectrum 
and the full distribution function of electrons, with arbitrary directional dependence. The latter was inte- 
grated out by Travis (1997) for isotropic electron distribution, which led to the following expression for the 
redistribution function: 

,_ "^l \ 1 _ e/(e» + e/)(27 + e»-e/) .^ 

2^i7^ \ yj\ - 2/isei/e/ -f- (e,/e/)2 2a,a/(aj + af) 






e^ef{l-fls) J (eje/(l-/is))2 V {a^af{a^ + aj)) 

2 2ai + af (oi + cif) — cuaf 2 Oj + 2a/ \ 



^'a^{a, + af) ^'''^^" (a^a/O^ ""^ajiui + af) 



(B3) 



Here, Oi = 1/(7 + ^i)^ - (1 + Ais)/(1 - fJ-s), af = y/il - &fY - (1 + Ms)/(1 - A*s), and r^ is the classical 
electron radius. In the Thomson limit ^ei <C 1. If one additionally assumes that Vi ^ Vj and 7 S> 1, the 
redistribution function reduces to 

w, = ^:L/i ^£ + ^1 \ (B4) 

which is equivalent to the result obtained by Reynolds (1982). 



B.l. Evaluating the Inverse Compton Emission Coefficient 

The multi-integral expression for the inverse Compton emission coefficient is frequently simplified using 
a (5-function approximation of the incident emission and/or the electron spectrum. We will demonstrate here 
that consideration of the broad-band incident spectrum is necessary. 

At a given location (i, x, y, z) the spectrum of the incident synchrotron emission in direction ili is given 
by Eqs. (28) and (29). Analysis of this solution shows that it can be represented roughly as 

^^t/^'''' '^*<^<'^^' (B5) 

exp(-j//j/<i), Vd < V, 

where vt, Vb, Vd are the self- absorption frequency, the break frequency, and the exponential drop-off frequency, 
respectively. We have tacitly assumed that ut < Vb < i^d, but one can readily visualize the shape of the 
resulting spectrum if the relations were different. The largest contribution to the incident synchrotron 
emission should be expected to arrive from directions that intersect the internal apparent position of the 
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shock front given by Eq. (11). The self-absorption frequency should increase only slightly due to an increase 
in the distance L between the internal apparent shock front position and the point [x, y, z) at which the 
emission coefficient is evaluated, vt oc h'^ly'^'^^) . The break frequency is expected to drop significantly, 
v\, ex i^^, until it reaches its minimum value when the apparent location of the shock reaches the geometric 
boundary of the excitation region. On the other hand, the drop-off frequency should remain constant and 
equal to v^ax = Cii37^jj^ as long as the apparent shock is inside the excitation region. Once the shock 
reaches the boundary the decrease in v^ will follow the changes in 7ma£c(Ai), where Ai is the time since the 
apparent shock crossed the geometric boundary. 

Although the redistribution function i]) has rather complex dependence on electron Lorentz factor and the 
incident and final frequencies, one can very crudely approximate it as ■0 ~ ^t I {'^^il'^) i^^ the Thomson limit. 
The innermost integral in the expression for the inverse Compton emission coefficient can now be estimated 
by substituting these approximations and performing the integration over Vi from vf /{2{1 — /is)7^) to vf. 
This gives 

{exp(-(7d/7)2), 7<7d 
f^r^S°, 7<i < 7 < 7fc ,T^«^ 

-{s-i)/2 _i (B6) 

v^ 7 , 7& <7<7t 



where the critical Lorentz factors (7^,76,7*) are found by solving 

"^ - 2(1 - /i,)72 



(B7) 



for 7 with V = Vd, Vb, and Vt- Function (B6) is then integrated over 7 from 7mi„(At) to ^maxi^t). From 
this solution it is apparent that electrons with Lorentz factors 7d < 7 < 76 and synchrotron seed photons 
with frequencies Vb < v < Vd provide equal contributions to the IC emission coefficient (see also McHardy et 
al. 1999). This is due to the fact that more abundant lower energy photons must scatter off of higher energy 
electrons, which are rarer, in order to contribute to the IC emission at Vf. Yet, less common higher energy 
seed photons can lead to the same final frequency Uf after being scattered by more plentiful lower energy 
electrons. 
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Fig. 1. — A cross-sectional diagram of a jet with a shock wave moving to the right and a stationary feature 
prior to colhsion. The moving shock is propelled by a fast stream of plasma in zone (1). Zone (2) contains 
the shocked quiescent jet material piled up in front of the fast flow. The shock between zones (2) and (3) 
moves at speed w^ in the AGN rest frame. The Mach disk slows down the quiescent flow [zone (3)], creating 
a dense region downstream [zone (4)]. The plasma that has been slowed down by the Mach disk is eventually 
reaccelerated to relativistic flow speeds in zone (5). The transition between (4) and (5) is probably gradual, 
but a well-defined boundary is used in the calculations for the sake of simplicity. 
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Fig. 2. — The space-time diagram of a collision between a stationary Mach disk and a shock moving down 
the jet with Lorentz factor F^ = 9; the speed of the quiescent flow in zone (3) is Fg = 4. Zones (2), (3), 
and (4) are the same as in Fig. 1 (the positions of the shocks shown in Fig. 1 corresponds to time before 
i' = 0). The collision occurs at time t' = at the location of the Mach disk. As a result of the collision a 
system of forward and reverse shocks is set up (solid lines for t' > 0). Two new zones, (r) and (/), contain 
the material originally from zones (2) and (4) that has been shocked by the reverse and forward shocks, 
respectively. Zones (r) and (/) are separated by a contact discontinuity (dashed line). The size H' of the 
excitation regions (in the AGN rest frame) for both forward and reverse shocks is indicated. Under the 
assumption that the plasma is fully relativistic in all zones, the resulting shocked plasma speed is v' — 0.64c, 
and the reverse and forward shock speeds are v'^ ~ 0.37c and v'f — 0.86c, respectively. 
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Fig. 3. — Shocked plasma parameters as a function of the moving shock Lorentz factor T'^ for ps = 1, p3 = 1, 
and Fg = 4; the plasma is assumed to be fully relativistic. Approximate solution for the shocked plasma 
Lorentz factor (dotted line) is shown in top panel. Speeds of f = l/\/3 and v — 1/3 are indicated by dotted 
lines in the middle panel for comparison (speed of light c — 1). Gradients along the jet axis may change the 
ratio of densities in zones (r) and (f ) . 
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Fig. 4. — A snapshot of the electron excitation structure while both forward and reverse shocks (vertical solid 
lines) are inside the excitation zones defined by Hf and H^- The shock fronts move in opposite directions in 
the rest frame of the emitting plasma at speeds Vf = 0.45c and Vr ~ 0.495c. The surface z — corresponds to 
the contact discontinuity between the forward and reverse zones. The emitting electrons are located between 
the forward and reverse shock fronts; gray scale corresponds to the maximum energy of electron, which is 
highest at the shock fronts. 
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Fig. 5. — A snapshot of the apparent excitation structure of the source as seen by a distant observer in the 
plasma rest frame; the line of sight subtends an angle Bobs to the z-axis. 
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Fig. 6. — Internal apparent shock orientation (solid lines) seen at an arbitrarily chosen point P in the forward 
zone at some time after the shock front passes it. The distance between point P and the apparent shock 
position at an arbitrary angle 9 to the 2;-axis depends on the light travel time from this position to point P 
and the speed of the shock front. 
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Fig. 7.-- The evolution of the synchrotron spectrum for 6obs = 0° in the rest frame of the shocked plasma in 
the forward shock region at times t/tac- 0.02, 0.2 (dot-dash), 0.5 (long dash), 1.0 (solid), 1.1, 1.5 (short dash). 
The bottom panel indicates the evolution of the synchrotron spectral index a^, defined by Fij{tobs) fx v"" . 
The following input parameters were used: R = H = 4 x 10^ s, v = 0.34c, B — 0.4 G, n = 10 cm^'', jmin = 
100, jmax = 2 X 10^, s = 2, Dq = lO'' Mpc. The synchrotron self-absorption frequency is z-'t = 8 x 10^ Hz 
and the maximum frequency is Vmax = 5 x 10^^ Hz. 
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Synchrotron Spectrum, 90° 
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Fig. 8. — The same as Fig. 7 but for 6obs — 90°, which corresponds to observation angle 6**^^ w 1/^1- The 
synchrotron self-absorption frequency is i/t — 10^° Hz and the maximum frequency is Vmax ~ 5 x 10^'* Hz. 
More complex behavior of the spectrum is apparent because the shock wave traces a more complicated 
geometrical pattern for this viewing angle. An additional spectrum at time t = 0.94 t^c is shown with a 
triple-dot-dashed curve. 
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Fig. 9. — Normalized synchrotron light curves for viewing angle 9obs — 0° at frequencies of 10^^ Hz (dashed), 
10^^ Hz (solid), and 10^'^ Hz (dot-dashed) in the rest frame of the shocked plasma. The time dependence of 
the spectral index for each light curve is shown in the lower panel. The ratio of the synchrotron decay time 
and the apparent crossing time is 2.0 (10^^ Hz), 0.7 (10^^ Hz), and 0.2 (10^^ Hz). The same input parameters 
were used as in Fig. 7. 
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Fig. 10. — The same as Fig. 9 but for viewing angle 9obs = 90°. The Ught curve profiles are roughly symmetric 
over a broad range of frequencies. The peak of the light curves at lower frequencies is delayed with respect 
to those at higher frequencies. The ratio of the synchrotron decay time and the apparent crossing time is 
0.8 (10" Hz), 0.3 (1012 Hz), and 0.08 (lO^^Hz). 
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Fig. 11. — The evolution of the inverse Compton spectrum for 9obs = 0° in the rest frame of the shocked 
plasma in the forward shock region. Compare with evolution of synchrotron spectrum in Fig. 7. 
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Fig. 12. — The same as Fig. 7 but for Qohs ~ 90°. Compare with Fig. 8. 
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Fig. 13. — Normalized inverse Compton light curves for viewing angle 6obs = 0° at frequencies of 10^^ Hz 
(dot-dashed), 10^* Hz (solid), and 10^" Hz (dashed) in the rest frame of the shocked plasma. The time 
dependence of the spectral index for each light curve is shown in the lower panel. 
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Fig. 14. — The same as Fig. 13 but for viewing angle 6obs = 90° 
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Fig. 15. — Forward (dotted) and reverse (dashed) synchrotron (10^^ Hz) and inverse Compton (10^^ Hz) 
light curves for viewing angle 9obs = 10° in the plasma rest frame. The time is normalized by tac,r = 3iac./- 
Solid curve is sum of both forward and reverse light curves. Dot-dashed curve corresponds to SSC emission 
from reverse region when the contribution of the seed photons from the forward region is neglected. We used 
the following parameters for the reverse region: H = 0.25 R, B = IG, n = 10 cm^'^; for the forward region: 
H = 0.05 i?, B = AG, n — 40cm~^. The other input parameters are the same for both regions: v = 0.34c, 

Imin = 10, -fmax = 2 X 10^, s ^ 2, Dq ^ 10^ MpC. 
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Fig. 16. — A model fit to the March 1999 flare in 3C 273; the observed flux density at 2keV and 2.2 /zm 
is indicated with filled circles. The observed and simulated X-ray spectral index variability is shown in the 
middle panel. Contributions from the forward (dotted line) and reverse (dashed line) regions are shown. 
A constant quiescent synchrotron flux of F^ ~ 82 mjy is used based on the average K-band flux detected 
before and after the flare. For the X-ray emission we use F^' = 10.3 /xJy and spectral index a^' — 0.67 to 
describe a constant baseline during the flare. We used the following input parameters: Hr = 0.6 x 10^ sec, 



Br 



0.7 G, Hr = 2.9 X 10^ cm-^ 7„ 



2 X lO'^ for the reverse region, and Hf ~ 0.15 x 10° sec, Bj — 2.0 G, 



nj = 4 X 10"^ cni^"^, jmaxj' = 4 X 10'^ for the forward region. The other input parameters were identical for 
both regions: R — 0.6 x 10^ sec, 7mm = 50, and s = 2.3. Modern cosmological parameters were used to 
convert the results of the simulation into the observer frame: Hq = 70 km s^^/Mpc, n.\ = 0.7. 
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Fig. 17. — Apparent crossing time tac for H = 2R (solid) and H = 0.2R (dashed) as a function of viewing 
angle in the plasma rest frame. We adopt R = I and v = 1/3. Dotted lines at tac = 7.4 and tac = 1-6 
correspond to the internal shock front crossing times tint for the central point (sec text). 
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Fig. 18. — Geometry for determining the apparent position of the shoek front. 



